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ABSTRACT 

We have obtained maps of the large scale outflow associated with the UCHII 
region G5.89-0.39 in CO and 13 CO (J=3-2), SiO (J=8-7,J=5-4), S0 2 (J=13 2 ,i2- 
13i,i3) an d H 13 CO + ( J=4-3). From these maps we have been able to determine the 
mass (3.3 M ), momentum (96 M km s -1 ), energy (3.5 x 10 46 erg), mechanical 
luminosity (141 L ), and mass loss rate (~lx 10~ 3 M Q yr -1 ) in the large scale 
outflow. The observationally derived parameters were used to guide 3D magne- 
tohydro dynamic models of the jet entrained outflow. Through the combination 
of observations and simulations, we suggest that the large scale outflow may be 
inclined by approximately 45° to the line of sight, and that the jet entraining the 
observed molecular outflow may have been active for as little as 1000 years, half 
the kinematic age of the outflow. 
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1. Introduction 

Over the past 20 years there have been significant advances in the study of star forma- 
tion. Several phases of the star formation process have been identified and characterized, 
beginning with a centrally concentrated core of molecular gas, which collapses to form a star 
surrounded by a proto-planetary disk (see for example, Shu et al. 1987). A large part of this 
progress has focused on the formation of low-mass stars, since they can form in isolation, are 
closer, and more abundant, reducing source confusion. Similarly, the formation timescales 
are long enough (tf or m > 5 x 10 5 yr, e.g. Hartmann 2000) for a considerable number of these 
objects to be detectable. 

Several studies have shown that low mass protostars have evolutionary phases both 
with and without outflows (Class 0/1 and Class II/III, e.g. Lada 1987, Andre et al. 1993) 
and that younger outflows appear to be more collimated and energetic (e.g. Bachiller & 
Tafalla 1999). A significant amount of effort has also gone into characterizing high mass 
protostellar phases but an analogous sequence has yet to be found. Meuller et al. (2002), 
and Shirley et al. (2002) used the same methods to determine physical properties in high and 
low mass star forming regions respectively and found accretion rates for high mass (M > 8 
M ) protostars were at least three orders of magnitude greater than those for low mass 
protostars. While in low mass systems, the Kelvin- Helmholtz timescale (t K H = GM^/R^L*) 
is too long to be of interest, in high mass systems it can become shorter than the free fall 
timescale, possibly causing the protostar to begin to radiate before fully accreting all of its 
material (Cesaroni 2004). Because high mass protostars are further away on average than 
their low mass counterparts, and are deeply embedded objects, they are generally studied 
through their effects on their environments, e.g., their molecular outflows. Studies of the 
outflow mechanisms can provide valuable insight into the accretion processes for low mass 
protostars since the two processes appear related (e.g. Andre et al. 1993). This could be 
similarly true for high mass protostars but the process may be complicated by such things 
as the radiation pressure and ionizing UV flux from the central star. 

The larger average distances, shorter accretion times and the clustered nature of massive 
star formation (e.g. Lada, Bally & Stark 1991) reduce our ability to isolate individual, 
massive, pre-protostellar cores. As a result, we still do not understand how massive stars 
form. Is it simply a scaled-up version of low-mass star formation in which accretion rates 
are high enough to overcome radiation pressure (e.g. McKee & Tan 2003; Cesaroni 2004), 
or do massive stars form from the coalescence of lower mass protostars (e.g. Bonnell et al. 
2001)? To investigate the outflow phenomenon in massive star formation, we will focus on 
the molecular outflow in G5. 89-0. 39 (hereafter G5.89) - also known as W28A2. 

The ultracompact HII (UCHII) region in G5.89 (cataloged by Wood & Churchwell 1989) 
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has a radius of 0.01 pc, is powered by an 05 ZAMS star (Feldt et al. 1999) and is expanding 
at a rate of ~ 35 km s^ 1 (at an assumed distance of 2 kpc, Acord et al. 1998). The 
molecular outflow in G5.89 was first studied by Harvey & Forveille (1988, hereafter HF88) 
who mapped the innermost square arcminute of the region and found that, at the edge of 
their map, high velocity line wings (Afpwzp ~ 25 km s _1 in 13 CO) were still apparent, 
suggesting the molecular outflow was much more extended than the area covered by their 
map. Recent sub-millimeter studies of this region have placed additional constraints on 
properties of the flow, such as momentum, kinetic energy, mechanical luminosity, age, mass 
loss rate, and force in both the ambient material and outflow lobes using the shock tracer 
SiO (Sollins et al. 2004; Acord et al. 1997, hereafter AWC97), showing this source to be 
quite powerful. The broad line wings, and extremely high mass powering source make this 
an interesting outflow to characterize. In terms of outflow energetics, the G5.89 outflow is 
the sixth most energetic in the Wu et al. (2004) study of high velocity molecular outflows, 
placing it in the top 2% of their outflows with energy calculations. 

In Section 2 we present the first observations of the full extent of the molecular (CO) 
outflow in G5.89, as well as SiO and serendipitous SO2 and H 13 CO + observations. We use 
these data to derive a number of physical properties of the outflow, which we use to constrain 
simulations of a jet entrained molecular outflow in Section 3. These high resolution mag- 
netohydrodynamic (MHD) simulations test the number of protostellar sources required to 
power the observed outflow, the inclination angle of the outflow, and the jet lifetime required 
to achieve the observed dynamics. Comparisons between observations and simulations will 
be drawn in Section 4. We then summarize our results in Section 5. For all calculations in 
this paper, we have adopted the more recent distance estimate to G5.89 of 2 kpc (i.e. Acord 
et al. 1998, Feldt et al. 2003), noting that this will cause systematic differences in derived 
properties from previous studies. 



2. Observations and Results 

Observations of the rotational transitions of 12 CO J=3-2, 13 CO J=3-2, SiO J=8-7 and 
SiO J=5-4 were taken in 2003 April and May at the James Clerk Maxwell Telescope (JCMT) 1 
on Mauna Kea, Hawaii. CO is a tracer of the bulk of the gas in the interstellar medium, 
since it is the second most abundant molecule, and the most abundant molecule with dipole 



1 The JCMT is operated by the Joint Astronomy Center on behalf of the Particle Physics and Astronomy 
Research Council of the United Kingdom, the Netherlands Organization for Scientific Research, and the 
National Research Council of Canada. 
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rotational transitions (e.g. Bachiller, 1996). We discuss our CO observations below in Section 
2.1. SiO is generally used to trace recently shocked gas, since silicon is generally frozen out 
of the gas phase in the ISM (e.g. Martin-Pintado et al. 1992). The transitions of SO2 
(J=132,i2-13i 5 i3) and H 13 CO + (J=4-3) were serendipitously observed during the SiO J=5-4 
and J=8-7 observations respectively, as discussed in Section 2.2. 

The central position of each map was the location of the Cesaroni et al. (1988) water 
maser at W28A2 (1); a B1950 = 17 A 57 m 26. s 803, 5 B1950 = -24°03'54."02. Feldt et al. (2003) 
recently identified a protostar at a 2 ooo = 18 h 00 m 30. s 44 ± 0. s 013, <5 200 o = -24°04'0".9 ± 0".2 
as the source of the UCHII region, but do not comment on whether this is the source of the 
large scale outflow. This is offset from the water maser by 2". 4, and does not correspond to 
any other water maser in the Cesaroni et al. list. Our map center lies south of the Feldt et 
al. (2003) center, in a region they show to be highly obscured by dust. All observations were 
reduced using the SPECX (Prestage et al., 2000) and CLASS (Buisson et al. 2002) software 
packages. 

Observations with Receiver A (230 GHz) were made in double sideband mode for SiO 
(J=5-4), had a half power beamwidth (HPBW) of 20", and a main beam efficiency (r] mb ) 
of 0.62. Observations made with Receiver B (345 GHz) were in single sideband mode for 
the 12 CO, 13 CO, and SiO J=8-7 observations, had a HPBW of 14", and r) mb = 0.62. All 
observations were taken in raster mapping mode, with the DAS configured to 760 MHz. 
Observations were coadded, first order baselines were subtracted, and the lines were fit with 
Gaussian line profiles. Table 1 provides a summary of our JCMT observations of the G5.89 
region. 

2.1. 12 CO and 13 CO 

We removed the Gaussian fit line centers from our spectra, leaving only the residual 
outflow component of the emission, the so-called line wing emission. Unless stated otherwise, 
when describing integrated intensities, we are referring to the integrated intensity of the 
residual line wing emission (either red shifted or blue shifted), not the emission from the 
line center. Gaussian models were fit to each 12 CO and 13 CO spectrum individually, using 
the Gaussian model routine within CLASS, and removed using the RESIDUAL command 
as shown for 13 CO at our map center in Figure 1. For reference, the average 12 CO Gaussian 
parameters with standard deviations were: Vlsr — 8.3±1.0 km s _1 , full width half maximum 
(FWHM) = 4.6±1.4 km s" 1 and = 26.7±11.1 K. The velocity extent of the wings was 
set as the range from -66 km s" 1 to 78 km s _1 , which is the full width at zero power towards 
the central position. Using this residual method, we have been able to separate the cloud 
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emission from the outflow emission. The 13 CO emission is very well fit by a Gaussian with 
line wings, suggesting that this would also be the case for the 12 CO if the self-absorption 
were not present. The blue side of the 12 CO emission does appear to be well fit with a 
Gaussian and line wing emission. Removal of the Gaussian from the 12 CO emission did 
tend to give negative intensities within the lower velocity red shifted absorption trough. To 
avoid incorporating these negative intensities into our subsequent calculations, we began 
calculating integrated intensities redward of 20 km s _1 . This velocity was chosen because it 
is well outside of the 1 a error bars on the average FWHM of the Gaussian, yet incorporates 
as much red shifted emission as possible. The presence of the absorption features will, 
unfortunately, tend to underestimate the integrated intensity of the line wing. 

To constrain better the properties of the G5.89 outflow over previous maps, our fully 
sampled 12 CO (J=3-2) map extends 3' along the flow axis, and 2' perpendicular to it, while 
our 13 CO (J=3-2) map covers the innermost 98" x 98" of the region. Figure 2 shows the 
total integrated intensity and the residual line wing emission of the 12 CO and 13 CO J=3-2 
towards G5.89. We define the edges of the outflow to be the locations in the CO map where 
the line wing integrated intensity is 10% of its maximum value. The maximum values are 
J Tdvb = 642 K km s -1 , J Tdv r = 452 K km s _1 for the blue and red wings respectively, 
while the rms noise is 6 K km s" 1 . From this definition, the 12 CO outflow subtends just 
under 2' (49" along the red lobe, and 56" along the blue) along the flow axis, and just under 
1' perpendicular to it and our maps encompass the entire outflow. In the less optically thick 
tracer ( 13 CO), the full extent of the outflow is 80" along the flow axis, and 50" perpendicular 
to it. At a distance of 2 kpc, the 12 CO outflow appears to extend 1.2 pc on the sky (without 
consideration for inclination angle effects). A literature search revealed no constraints on 
the inclination angle of this outflow. The outflow does however appear extended on the sky 
with little overlap between the red and blue lobes, thus it is likely that the outflow is not 
primarily oriented along the line of sight (i.e. 90° inclination). We compare our observations 
to the models of Cabrit & Bertout (1990) in Section 2.3. 

Every position in the 12 CO map contains absorption due to cold clouds along the line 
of sight at Vlsr=13.7 km s _1 and 20.2 km s" 1 . The presence of these absorption features 
affects the total observed integrated emission from the red shifted outflow which produces 
uncertainties in calculations requiring integrated intensities. Figure 3 shows spectra from 
the central position in each tracer and the two line of sight absorption features in 12 CO are 
clearly distinguishable. Figure 1 of HF88 shows three CO isotopes (CO, 13 CO, and C 18 
in the J=l-0 transition) and two CS isotopes (CS, and C 34 S in the J=2-l transition). In 
both HF88 and our study, 12 CO is the only molecular species which appears to suffer from 
the two absorption features. Thus we assume our 13 CO observations are free of line-of- 
sight absorption. Note, however, that while Figure 2 shows a red outflow lobe which is less 
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extended than the blue lobe, this is not likely to be an effect of the absorption features. 
Since the red absorption features appear at every position in the map, calculating the extent 
of the outflow as a percentage of the maximum causes the absorption effect to cancel out. 
The same is not true of calculations of physical quantities, however (see Section 2.3). 

Figure 4 shows the position- velocity (PV) diagram of 12 CO emission along the flow axis 
(AS = 0). It reveals that the highest velocity gas is concentrated towards the center of the 
outflow, suggesting that the gas decelerates away from the source (see for instance Cabrit & 
Bertout 1986, 1990). The solid lines show a v oc r _1 trend, which the bulk of the gas tends 
to follow. For reference, the dashed line shows a v oc r trend, which represents accelerating 
gas. At very low levels on the blue shifted side of the PV diagram, there does appear to be 
some accelerating gas in this region. We suggest that this may be a contribution from the 
accelerating outflow in this region detected by Sollins et al. (2004). Their outflow (in SiO) 
appears to have a velocity extent of ~ 50 km s -1 and an angular extent of 15". Our tenuously 
detected second (accelerating) outflow appears to extend 60" in CO, beyond their primary 
beam. However, since we do not have the spatial resolution to confirm this hypothesis, we 
do not discuss it further. 



2.2. SiO, S0 2 and H 13 CO+ 

Silicon is a key constituent of dust grains, and the passage of a shock wave can remove 
Si from dust grains via sputtering (e.g. see Field et al. 1997). The elevated temperatures 
arising from a shock can also liberate Si-bearing species from frozen grain mantles. Once in 
the gas phase, Si atoms react quickly with oxygen, forming SiO. SiO, however, can also be 
quickly oxidized into Si02, or freeze back onto grain mantles in < 10 4 years (e.g. Pineau des 
Forets et al. 1997). Thus, SiO is an excellent tracer of recent shocks, such as those produced 
by molecular outflows (e.g. Bachiller 1996). The combination of CO and SiO observations, 
therefore, can yield a broader understanding of the G5.89 outflow. 

Figure 5 shows 80" x 80" maps of the SiO J=5-4 and J=8-7. In the SiO J=5-4 data we 
suggest that the serendipitous line seen at Vlsr ~ 77 km s^ 1 is SO2 J=132,i2-13i,i3 in the 
upper sideband (v = 225.154 GHz), as seen in a 215-247 GHz line survey of Orion A (Sutton 
et al. 1985). Models of shock propagation and chemical evolution concerning sulfur bearing 
molecules show that at temperatures in the 100-300 K range, S0 2 is the dominant sulfur 
bearing species (Doty et al. 2002, Charnley, 1997). The S0 2 emission in our observations 
appears to have the same spatial extent as the SiO emission, suggesting a common origin 
(shocked, warm gas). Another serendipitous line was detected in the SiO J=8-7 observations, 
and is shown in the inset of Figure 3. This line is H 13 CO + J=4-3 at v — 346.999 GHz, 
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consistent with the molecular line survey of G5.89 by Thompson & MacDonald (1999). 

2.3. Physical Parameters Derived from Observations 

Assuming that the 12 CO emission is optically thick, that the level populations are in 
LTE, and a beam filling factor of one, we derive a kinetic temperature at the central position 
of ~ 80 K from the peak brightness of the line. We note however, that one transition of 
CO is not a sensitive temperature tracer (especially of material deeper into the cloud than 
the r ~ 1 surface) and suggest that this temperature is a crude lower limit to the actual 
kinetic temperature of the region. Observations of NH 3 were used by AWC97 to place a 
lower limit on the temperature of the SiO emitting region at 100 K. The presence of SO2 in 
our observations also suggests a lower limit of Tk ~ 100 K (Doty et al. 2002). Thus, we 
have adopted a kinetic temperature of 100 K for the central position of G5.89. 

Schilke et al. (1997) presented a set of C-type chemical shock models with a variety 
of different ambient densities and shock velocities. Their results include level populations 
up to SiO J up = 15 (see their Figure 6) as well as integrated intensities and line ratios for 
three of the SiO rotational transitions observable by ground based telescopes (J=2-l, 5-4, 
and 8-7). Convolving our J=8-7 observations to the beamsize of our J=5-4 observations, we 
find a line integrated intensity ratio of [8-7]/ [5-4] = 0.83. Comparing this to the Schilke et 
al. model, we estimate n H = 10 7 cm~ 3 and t> s =28 km s -1 , but note that at v s > 30 km s _1 , 
this ratio is only weakly dependent on shock speed. A simple large velocity gradient (LVG) 
model of the SiO emission from the two observed transitions, however, gives similar ambient 
densities (n = 1.4 x 10 7 cm -3 ). These densities are higher than those derived by previous 
studies (n « 10 6 cm -3 ; Plume et al. 1997, AWC97). 

The observed outflow lobes in G5.89 are not highly inclined towards the line of sight, 
and a comparison of Figures 2 and 4 to the models of decelerating outflows of Cabrit & 
Bertout (1990), show that this outflow lies somewhere between their i = 10° and i = 50° 
cases, closer to the i = 50° case. We suggest an inclination angle of approximately 45° with 
respect to the plane of the sky. 

The dynamical properties of the outflow (such as momentum, energy, luminosity, and 
age) are best studied when the line emission is broken down into velocity intervals, allowing 
for examination of the mass and dynamics in each velocity bin. Summing over all the bins 
gives the total momentum, energy, etc., in each outflow lobe. Figure 6 shows channel maps 
of the 12 CO emission from G5.89, while Table 2 shows the dynamical properties calculated, 
the equations used, and values obtained for each outflow lobe. The kinematic age of the 
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outflow was determined by dividing the extent of the outflow (scaled by a factor of sin 45° 
for the inclination angle) by the average gas velocity (scaled by a factor of cos 45° for the 
inclination angle), resulting in an age of tkm — t^° 3 tan 45° ~ 2000 yr. This age is consistent 
with that given by Wu et al. (2004), who used similar methods to derive outflow lifetimes. 

Using the column density estimated from the 13 CO observations, an assumed 12 CO/ 13 CO 
abundance ratio of 55 (Langer et al. 1984), and the derived abundance of H 2 in this region 
with respect to CO (1.7xl0 4 ; van der Tak et al. 2000), we can determine the abundances 
of our other molecular species. Based on our SiO J=8-7 observations, which have a higher 
signal to noise ratio than our J=5-4 observations, and the same LTE approximation used 
to determine the 13 CO column density (such as an ambient temperature of 100 K), we find 
[SiO]/[H 2 ] = 3xl0 -10 . This abundance is consistent with SiO being recently released from 
grain mantles into the gas phase (Schilke et al. 1997). Using the RADEX online calculator 2 , 
we were able to determine the column density of our S0 2 gas (N S o 2 ~ 3 x 10 15 cm -2 ), 
and an abundance ratio of [S02]/[H 2 ] = 2xl0~ 8 . We compare this S0 2 abundance to the 
hot core sulfur models of Charnley (1997), to find that S0 2 reaches this abundance after 
approximately 1900-2500 yr. This age is consistent with the kinematic age of the outflow 
described above. 

If we assume an accretion rate 10 times the (constant) outflow rate (shown in Table 
2), it would take 2000 years for 60M Q (the mass of an 05 star) to accumulate. Studies of 
accretion and outflow rate ratios suggest that even higher ratios exist (e.g. Beuther et al. 
2002), and thus it is possible for the accretion timescale to be shorter than 2000 yr. 

We can easily compare our results to those of Beuther et al. (2002). They determined 
physical parameters for 33 massive molecular outflows using methods similar to ours. They 
did use a slightly higher 13 CO abundance (8.9xl0 5 ) than we did (3.2 xlO 5 ), but when we 
correct for this, we find that our mass, momentum and mass loss rates are lower than their 
average values, but that our energy, luminosity and force are all significantly higher than their 
average values. There are only two sources in their source list (18264-1152 and 19410+2336) 
which have higher outflow energies than G5.89, and only one (19410+2336) which has a 
higher mechanical luminosity and force. Our derived values for this source are consistently 
smaller than those derived by previous studies (i.e. HF88, AWC97, Sollins et al. 2004). 
The authors of those studies each used the total integrated intensity of their observed lines, 
whereas we only used the outflowing line wing emission to determine outflow properties. 

To further the analysis of the G5.89 flow, we used the observationally determined proper- 
ties of the outflow to guide magnetohydrodynamic (MHD) simulations of a disk-wind driven 



2 http://www.strw.leidenuniv.nl/^moldata/radex.php 
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molecular jet which entrains the ambient molecular gas into an outflow. This allows us to 
produce an outflow mimicking our observations, from which we can draw more information, 
and conclusions about the nature of massive star formation. 



There are a number of theories suggesting entrainment mechanisms of molecular out- 
flows associated with massive star formation. While we do not discuss the specifics of the 
driving mechanism, we assume an accretion model, and that a disk wind driven molecular 
jet (e.g. Pudritz & Norman, 1986, Raga & Cabrit 1993, Masson & Chernin 1993) entrains 
the surrounding molecular gas. Regardless of the specifics of jet generation and outflow 
entrainment, we can estimate the inclination angle, kinematic age, and length of time the 
jet is powering the outflow by comparing our MHD simulations to observations. We can use 
observational constraints as initial conditions for the simulations. 

There are two distinct sets of input parameters required for our simulations: source and 
ambient. The source parameters were based on an 05 zero age main sequence (ZAMS) star, 
like the one powering the G5.89 outflow (Feldt et al. 1999). Ambient parameters were fixed 
based on the observational analysis presented above. All of these parameters are given in 
Table 3. Parameters derived from these quantities are discussed below. 

Within our model, the outflow particles are attached to the magnetic field lines, and 
so, are not launched from the surface of the protostar, but from the magnetic footpoints 
in the disk (r D ). They are then accelerated until they reach the Alfvenic radius (r a ) where 
collimation of the outflow occurs. It is the ratio between these two radii (r /r a ) which 
determines the amount of collimation and acceleration within the jet (e.g. Blandford & 
Payne 1982, Pudritz & Norman 1986). We will define this ratio to be a 2 , resulting in a 
becoming the scaling factor between the Keplerian velocity at the edge of the star, and the 
point at which the particles are launched into the outflow (Ouyed & Pudritz, 1997): 



Setting a = 0.4 gives v w = 552 km s" 1 . This velocity is only slightly slower than the Ha 
line wings in HH 444 discussed in Andrews et al. (2004). As a control, we also set a = 0.2, 
to see the effects of lower wind velocities on the final extent of the outflow. 

In general, how much of the source bolometric luminosity is transfered into the outflow 
mechanical luminosity is a poorly constrained parameter. There are a number of ways of 



3. 



Simulations and Results 
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dissipating the energy produced by a protostar, such as momentum transfer to the outflow, 
shock heating of the ambient gas and atomic and molecular line cooling. Since the source is 
likely an 05 ZAMS star, we can constrain the bolometric luminosity to be L\ M \ = 8 x 10 5 
L (Carroll & Ostlie, 1996). Similarly, from our outflow observations, we can constrain the 
mechanical luminosity of the flow to be -L mcc h = 140 L Q (see Table 2), and can express it as a 
fraction of the source bolometric luminosity (L mech = /3L bol ). Other authors have suggested 
the possibility of multiple sources powering the outflow in G5.89 (i.e. Feldt et al. 2003, 
Sollins et al. 2004), however with our (3 10~ 4 we find, as do they, that we do not require 
a second source to account for the energetics of the outflow. 

We observationally constrained (see Section 2.3) the ambient density within the central 
14" of the G5.89 outflow to be ~ 10 7 cm -3 , and applied a Plummer profile covering a region 
of radius 0.6 pc to mimic density gradients common to protostellar regions, which results in 
an average density of 10 4 cur 3 (i.e. Whitworth & Ward-Thompson 2001, Boily & Kroupa 
2003). To simplify our initial simulations, we used the averaged density (10 4 cm -3 ) for the 
entire region, with ambient density gradients being left to future study. A summary of all of 
our observationally constrained initial conditions is shown in Table 3. 

Further simplifying assumptions include use of only one molecular line coolant, CO, and 
that it is not destroyed by the jet bow shock. CO is a coolant inherent to the code, and at 
the temperatures suggested by observations, is one of the strongest coolants in the cloud (e.g 
Smith & Rosen 2003). Simplyfing the cooling may slightly inflate the CO line profiles, but 
not significantly. We also do not treat the heating, destruction and reformation of CO by the 
jet as we are more interested in the extent and lifetime of the jet and surrounding outflow, 
than the chemistry produced by the jet. The momentum transfered to the ambient medium 
is the quantity we are attempting to measure, which is independent of the nature of the 
medium (e.g. whether it is atomic or molecular). Neither of these simiplifying assumptions 
should adversely effect our results. 



3.1. Computational Details 

Simulations were conducted using the Zeus-MP astrophysical fluid dynamic code (Nor- 
man 2000) on the CAPCA 3 computer cluster. The cluster consists of 64, 2.4 GHz, Linux 
based processors connected by a 1 GB network. Results of these simulations were imaged 
and analyzed using JETGET (Staff et al. 2004), an interactive data language (IDL) based 



3 Animations of these simulations (like the panels shown in Figure 7) can be viewed by following the 
"Animations" link at www.capca.ucalgary.ca. 
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visualization program for use with hierarchical data format (HDF) files produced by Zeus 
like codes in two or three dimensions. 

The simulations were run in three Cartesian coordinates, with the protostar comprising 
the central pixel (or grid zone) of the simulated region. The grid used for these simulations 
was 210 x 240 x 210 pixels, with a resolution of x\ x x 2 x x 3 = (1.1 x 5 x 1.1) x 10~ 3 pc, 
giving a total extent of 0.24 x 1.2 x 0.24 pc in the plane of the sky. For reference, the 
X\ plane is perpendicular to the outflow axis, and corresponds to the line-of-sight in the 
observations. The x 2 plane is parallel to the outflow axis and corresponds to right ascension 
(in projection), while the x% plane is perpendicular to the outflow axis, corresponding to 
declination (in projection) for the large scale G5.89 outflow. 

3.2. Simulation Results 

To determine how long the jet actively powers the outflow, we ran a number of simula- 
tions with varying jet activity timescales. For each simulation, the jet was turned off at t g , 
and the molecular outflow was allowed to evolve to the kinematic age of the system (t^m = 
2000 yr; Section 2.2). Simulations were run using t s = 300 yr, 500 yr, 1000 yr, and oo (jet 
does not shut off within i kin ). We then compared the size of the simulated outflow lobes at 
i kin = 2000 yr to the size of the observed CO lobes (1.2 pc at a distance of 2 kpc uncorrected 
for inclination angle; Section 2.1). The shorter jet lifetimes (t oS = 300 and 500 yr) were 
used to test the simulations and ensure we were not getting jet activity timescales which 
were unreasonably short. Given the ambient density of 10 4 cm" 3 , the outflows powered by 
short lived jets (t Q g = 300 and 500 yr) reach full extents of 0.7 and 0.9 pc, respectively. As 
expected, the spatial extents, and gas velocities of these two sets of simulations were too low 
when compared to the observed outflow in G5.89. 

We found that setting a = 0.2 did not produce results consistent with the G5.89 outflow. 
With a = 0.2, an outflow with a jet lifetime of t a g = oo could not reach the extent of the 
observed outflow. Thus, a was set at 0.4 for the rest of our simulations. We acknowledge 
the possibility that a ^ 0.4. However, for a < 0.4, the simulated outflow will not reach the 
observed spatial extent of the outflow. Higher wind velocities (a > 0.4) would decrease the 
jet activity timescale. Observational evidence for higher velocity winds from protostars is 
scarce. 

Figure 7 shows the results of simulations with longer t Q s, at various epochs in the outflow 
evolution. The x 3 plane (y axis) has been expanded (with respect to the x 2 plane, or x axis) 
in order to show the details of the jet. The first five panels show the evolution of the t oS 
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= 1000 yr case in 400 yr steps, while the sixth panel (shown with a different color scheme) 
shows the final result for the t a s = oo case. Both sets of simulations appear to reach the 
same spatial extent as the observed outflow. For the longer t g- cases, it becomes harder to 
differentiate between simulation sets. In both cases, the outflow appears to have the same 
extent, with the highest velocity gas varying by less than 1% between simulations. The 
collimation factor of the outflow in both cases is £\ cn gth /^width ~ 3, consistent with other 
high mass protostellar outflows (see for example, Beuther et al. 2002). The highest velocity 
outflow gas in both cases is ~ 70 km s -1 , with the jet velocity peaking at ~ 430 km s -1 . 
A noticeable difference between the simulations is the extent of the underlying jet. In the 
t Q g = 1000 yr case, the latest time panel of Figure 7 (bottom, middle) only shows the jet 
towards the ends of the outflow. In the t Q R = oo case, the jet extends from the central source 
to the edges of the outflow. 

Figure 8 shows the density of the simulations both at the resolution of the simulations, 
and convolved to the resolution of the JCMT beam at 345 GHz (simulated pixels convolved 
with a 15" Gaussian). Unlike Figure 7, where the y axis was stretched to show the details 
of the jet, the x and y axes in this figure have the same scaling. The outflow lobes appear 
asymmetrical for two reasons: First, they represent single channel maps centered at 4 km 
s _1 (with a width of 4 km s _1 ) with respect to the source velocity. This means that the 
emission in the red outflow lobe dominates the emission in the blue outflow lobe. Second, we 
have introduced and inclination angle of 45° as described in Section 2. This inclination angle 
appears to match well with observations. It is clear from Figure 8 that higher resolution 
images are required in order to resolve the jet, both spatially and kinematically. 

Figure 9 shows the density of the jet at the same times shown in Figure 7 (with the first 
panel, t kin = 400 yr, removed), 0.2 pc from the source along the jet axis, one third of the 
way from the source to the edge of the simulation. We can see that as the jet propagates 
through the medium, the bow shock entrains the surroundings, clearing out a cavity around 
the jet. The density of the jet was allowed to vary as a free parameter based on the ambient 
density and the force imparted by the central source. The jet itself appears to have a density 
of 4xl0 4 cm -3 , with the swept up surrounding material having a much lower, but non-zero, 
density. As the simulation progresses, Figure 9 shows that the bow shock propagates further 
from the jet axis, causing the outflow to lose collimation (at that point) with time. 

4. Discussion 

To determine the inclination angle of the outflow in G5.89, we compared our observations 
to the models of Cabrit & Bertout (1990), as well as our own simulations. In both cases, we 
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found the inclination angle to be ~ 45°. We can apply this correction factor to our derived 
outflow properties (as was done in Table 2), as well as the total extent of the outflow. 
Correcting for the inclination angle changes the total extent of the observed outflow to be 
1.6 pc. Applying the same correction to the simulated outflow, the extent of the outflows 
in the t g = 1000 yr and oo cases is 1.55 pc, similar to that observed. The same is not 
true of the shorter jet lifetime simulations. Even with inclination angle effects taken into 
consideration, the extent of the outflows in the t a s = 300 yr and 500 yr cases (respectively) 
were 1 and 1.3 pc, too small to match observations. 

From our simulations, the only distinguishing characteristic between the t Q ff = 1000 yr 
and oo cases is the length of the jet, which we cannot observationally constrain. In the t g 
= 1000 yr case, the jet is only visible towards the edges of the flow, whereas in the t off = oo 
case, the jet is visible from the source to the edge. It is interesting to note that the jets 
with lifetimes of at least 1000 yr both appear to reach the same extent within 2000 yr, the 
kinematic age of the outflow. It is unclear why this happens, however we suggest it may be 
a result of the chosen ambient density. At lower densities, the jet would be able to easily 
entrain the surrounding gas, while maintaining a high enough velocity to reach larger sizes 
on much shorter timescales. Conversely, with higher ambient densities, the jet transfers its 
momentum to the outflow and cannot reach the extent of the observed G5.89 outflow. This 
issue will be addressed in a future paper. 

A comparison of our observations with our simulations suggests a jet lifetime of > 
1000 yr. This could indicate that the large scale molecular outflow in G5.89 is a remnant 
of previous star formation activity. If we use the accretion to outflow rate ratio discussed 
earlier (10; Hartmann, 2000), in order for 60 M to buildup requires 2000 yrs. However, if we 
increase the ratio to 14 (a rough upper limit suggested by Beuther et al. 2002), the outflow 
could conceivably be powered by a jet active for as little as 1500 yr. Since jet activity is 
linked to accretion, this would give an accretion timescale of only 1500 yr. 

Some theories for the formation of UCHII regions suggest that they can only form once 
accretion has fallen below a critical value (e.g. Yorke 1984, Garay & Lizano 1999, Yorke et 
al. 2002). The UCHII region surrounding G5.89 has been previously shown to be 600 yr 
old (Acord et al. 1998, Zijlstra & Pottasch 1988), and physically linked to the source of the 
outflow (Sollins et al. 2004). If these results are correct, and accretion must halt before the 
formation of the UCHII region, an accretion timescale of 1500 yr appears to fit with the age 
of the UCHII region and kinematic age of the outflow (2000 yr). 

Using Figure 8, we can compare the column density in the observed outflow to that in 
the simulations. The H2 column density within the central convolved pixel of the simulations 
is 2.7xl0 20 cm -2 integrated along the emitting region of the simulation and the 4 km s~ 1 
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velocity bin centered at 4 km s" 1 . While this portion of the observed red shifted CO is self 
absorbed, we can use the equivalent velocity bin on the blue shifted side, as the simulations 
should be symmetric. Using the integrated intensity of that CO bin (using the same method 
described in Section 2.3), the CO column density is ~ 1 x 10 17 cm ~ 2 . This results in a 
CO abundance of [CO]/[H_2] = 3.7xl0~ 4 , which is approximately within a factor of two of 
the derived abundance for this region from van der Tak et al.(2000), suggesting that the 
observations and simulations are well matched. 

Other signposts of ongoing star formation include the presence of a protostellar disk, 
and bullets of molecular gas near the central source detectable in the infrared. To date, no 
disk has been detected in G5.89 (Sridharan et al. 2002). If molecular bullets were released 
from the central source within the last 500 yr (t kin = 1500 yr), we should be able to detect 
them within a given radius of the source. If we assume they would be moving outward at 
the same rate as the highest velocity outflow (not jet) gas (70 km s _1 ), the farthest they 
could reach in 500 yr is 0.03 pc. Any molecular bullets within that radius would have been 
emitted more recently, and suggest that the jet has been active for longer than 1500 yr. The 
field of view of the infrared observations of Feldt et al. (2003) is also approximately 0.03 pc, 
and no bullets can be seen in their observations, however this region is highly extincted. 

Given the results of our observations and simulations, and lack of evidence to the con- 
trary, it is plausible to suggest that the outflow in G5.89 is what remains of a now extinct 
jet. The results suggested here have interesting implications for the study of massive star 
formation. It is possible that accretion timescales are much shorter than previously believed 
(as little as 1000 yr). It may also be possible that the apparent association between massive 
outflows and UCHII regions may be an artifact of the short massive star formation timescale; 
that an observed remnant outflow may not be indicative of ongoing accretion. While others 
have suggested that infall can continue through a Hypercompact HII region (e.g. Keto 2003), 
we suggest that this is not the case for G5.89. 

5. Conclusions 

We have created the largest map of the G5.89 outflow to date, and have, for the first time, 
captured the full extent of the large scale outflow. Through comparison of simulated and 
observed channel maps, we suggest the outflow is inclined by ~ 45° to the line of sight. This 
implies that the combined extent of the two outflow lobes is 1.6 pc. From our observations, 
we were also able to constrain the mass, momentum, energy, etc. of the outflowing gas. 

We used these observationally derived parameters as input constraints on a set of MHD 
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simulations of a jet entrained molecular outflow. Our results suggest only one source powered 
the outflow, and that the jet must have been active for a minimum of 1000 yr to account 
for the spatial extent of the observed outflow. This jet activity timescale is consistent with 
a 600 yr old UCHII region, and a 2000 year old outflow. 

The authors would like to thank T.K. Sridharan for helpful discussions, and the anony- 
mous referee for strengthening the paper. We would also like to acknowledge the National 
Science and Engineering Research Council of Canada for their financial support. 
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Tracer Map Parameters System Parameters 



Transition 


V 


Mapsize 


Spacing 


T 

± sys 


T 

± rms 


T 


Av 


tint 




(GHz) 


(Pixels) 


(") 


(K) 


(K) 




(km/s) 


(s) 


i2 C J=3-2 


345.796 


28x17 


7 


507.9 


0.241 


0.20 


1.08 


2.57 


13CO J=3-2 


330.588 


15x15 


7 


759.2 


0.190 


0.37 


1.13 


1.27 


SiO J=8-7 


347.331 


12x12 


7 


660.4 


0.059 


0.26 


1.08 


0.97 


SiO J=5-4 


217.105 


8x8 


10 


467.3 


0.190 


0.15 


3.45 


1.82 



Table 1: Observation Parameters. System parameters (system temperature, rms noise limit, 
atmospheric opacity, velocity resolution and integration times) are given for the central 
position in each map. 







Outflow Value 


Characteristic 


Equation 


Red blue 


Mass (M Q ) 


M = N H2 -A-m H jM Q 


0.6 2.7 


Momentum (M Q km s~ 4 ) 


P = Y.i{™,i/M Q )\v i -v LSR \ 


17 79 


Energy (xlO 45 erg) 


E = 1/2 -v LSR ) 2 


5.6 29 


Luminosity (L Q ) 


Lmech Eft 


23 118 


Mass loss (xlO" 3 M yr _1 ) 


M = M/t 


0.3 1.4 


Force (xlO~ 3 M Q km s _1 yr _1 ) 


P = P/t 


8.5 59 



Table 2: Flow characteristics corrected for an inclination angle of 45°. The t in the mass 
loss rate and outflow force equations refers to the kinematic age of the outflow derived from 
observations (tkin = 2000 yr). 



Simulation initial conditions 



Source Parameters 


Ambient Constraints 


M* 60 M 


T ism 100 K 


R* 12 Rq 


n ism 10 4 cm -3 


v w 552 km s _1 


c s 1 km s~ x 


P 10" 4 





Table 3: Constraints on values for source and ambient parameters. 
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Fig. 1. — Gaussian fit removal from the 13 CO spectrum at our map center, (a): Origi- 
nal 13 CO spectrum from our map. (b): Gaussian fit (dashed line) and 13 CO spectrum, 
(c): residual 13 CO emission once the Gaussian profile was removed using the RESIDUAL 
command in CLASS. 




Fig. 2. — 12 CO (left) and 13 CO (right) integrated intensity maps overlaid by red (dashed) and 
blue (solid) line wing emission. Contours start at 10% of the integrated line wing intensity 
and increase in 10% increments. The grey scale shows total integrated intensity in units of 
K Km s -1 . 
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Fig. 3. — Spectra from the central position in each observed 12 CO, 13 CO, and SiO transition. 
Note that the 13 CO J=3-2 and SiO J=5-4 and J=8-7 spectra are offset from the 12 CO by 
18 K, 25 K, and 35 K (respectively). The redshifted material in the 12 CO observations 
is absorbed at two different sets of velocities by cold line of sight clouds. Based on the 
observations of HF88, we suggest this is not the case for the 13 CO and SiO. The SiO spectra 
have been multiplied by a factor of 10 for comparison with the CO spectra. 
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Fig. 4. — P-V diagram for G5.89. The two solid white lines represent contours of v oc r _1 , 
while the dashed diagonal line represents v oc r. There appears to be evidence for both 
a larger scale decelerating outflow and smaller scale accelerating outflow. The slice shown 
above corresponds to A DEC = in Figure 2. 
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Fig. 5. — Total integrated intensity maps ( f T^dv) for each non CO transition observed at 
the JCMT. Specific molecular transitions are given in the upper left hand corner of each 
map. The maps are on the same spatial scale, but the integrated intensities represented 
by the grey scale are as shown at the top of each map (in units of K km s _1 ). For maps 
of integrated SiO J=5-4, S0 2 J=13 2 ,i2-13i,i 3 , SiO J=8-7, and H 13 CO+ J=4-3, the contours 
begin at 10, 5, 5, and 10 K km s" 1 respectively, and increase in steps of 10, 5, 10 and 10 K 
km s -1 respectively. 
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Fig. 6. — Channel maps of the G5.89 outflow in 12 CO J=3-2. The number in the upper right 
hand corner of each frame is the middle velocity of the bin (4 km s _1 per bin) with respect 
to the Vlsr of the source (9 km s _1 ). 
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Fig. 7. — The first five panels (with green backgrounds) show velocity contours, and the 
evolution of the simulated outflow in which the jet was shut off after 1000 years. Panels 
show 400 year intervals, starting at 400 years and ending at 2000 years. The bottom right 
panel (with red background) shows the t oS = 2000 yr outflow at a kinematic age of 2000 yrs 
for comparison to the t Q s = 1000 yr case (bottom center). The extent of the two outflows 
are the same, however, the jet can be traced back to the central source in the t fr = 2000 yr 
case. The two dimensional slices of the three dimensional simulations are in plane (or slice) 
105 out of 210 slices, through the central source. The halftone scale represents the velocity 
of the gas, with the minimum (~ -430 km s _1 ) and maximum (~ 430 km s" 1 ) velocity given 
for each panel. The scales in the X2 and x% planes are in units of 0.1 pc. 
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Fig. 8. — Column density contours for simulated jet and outflow, shown at an inclination 
angle of 45°. The maps shown are channel maps, taken over the simulated velocity range of 4 
km s _1 centered at 4 km s _1 (e.g. 4±2 km s _1 with respect to the source velocity). The top 
panel shows the jet, with equal scaling on the x and y axes, while the bottom panel shows 
the jet convolved with a 0.145 pc (15" at the distance of G5.89) Gaussian for comparison to 
our JCMT observations (e.g. the panel labelled as "+3" in Figure 6). There are ten column 
density contours in each panel, starting from 1.2 x 10 21 cm~ 2 and 2.7 x 10 20 cm -2 (for the top 
and bottom panels respectively) and decreasing in intervals of 10% of the maximum value. 
The H 2 column density in the central pixel of the lower plot is only a factor of two off from 
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Fig. 9. — Jet density is shown here as a function of time. Each panel represents 800 (top 
left), 1200 (top right), 1600 (bottom left) and 2000 (bottom right) years into the t ff = oo 
case, at a distance of 0.2 pc from the source along the jet axis. The density on the y axis 
is given in units of 10 4 cm -3 , while the units on the x axis are given in units of 0.1 pc from 
the jet axis. 



